Vibrational density of states of silicon nanoparticles 
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The vibrational density of states of silicon nanoparticles in the range from 2.3 to 10.3 nm is studied 
with the help of molecular-dynamics simulations. From these simulations the vibrational density of 
states and frequencies of bulk-like vibrational modes at high-symmetry points of the Brillouin-zone 
have been derived. The results show an increase of the density of states at low frequencies and a 
transfer of modes from the high-frequency end of the spectrum to the intermediate range. At the 
same time the peak of transverse optical modes is shifted to higher frequencies. These observations 
are in line with previous simulation studies of metallic nanoparticles and they provide an explanation 
for a previously observed discrepancy between experimental and theoretical data [C. Meier et ai, 
Physica E, 32, 155 (2006)]. 

PACS numbers: 63.22.Kn,62.25.Jk,81.05.Cy,02.70.Ns 
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I. INTRODUCTION 

Modern nanotechnology is leading to the development 
of a new generation of applications and devices based on 
nanostructured materials. An important part of the de- 
velopment of many nano-devices is the thermal design. If 
the heat generated inside a device cannot be transported 
away efficiently, the system will fail or even break due 
to overheating. Thermal properties like heat capacity 
and thermal conductivity as well as many other material 
properties are strongly influenced by the vibrational den- 
sity of states (VDOS). For this reason a profound under- 
standing of the laws governing the vibrational properties 
of nanostructured materials is of high technological and 
fundamental interest. 

This article focuses on the vibrational properties of 
silicon nanoparticles. Silicon is an important engineering 
material whose bulk properties have been studied exten- 
sively. Its favourable electrical properties have made it 
the base material for all microelectronics and advanced 
technology exists for the manufacturing of nanostruc- 
tured silicon based devices. 

The vibrational properties of bulk crystalline Si have 
been studied both experimentally using neutron scatter- 
ing techniques^ as well as theoretically using empiri- 
cal and first-principles methods.^— In addition to this, 
Raman scattering experiments have been used in nu- 
merous studies of the vibrational properties of silicon 
nanostructureSf&~- e In many of these studies it is found 
that the first order peak due to the transverse optical 
(TO) phonon at the T point is broadened and its maxi- 
mum shifted to lower energies. The phonon-confinement 
model&ii attributes this behaviour to the fact that in 
nanostructures more vibrational modes can become Ra- 
man active than in bulk crystals since the translational 
symmetry is broken and consequently no Ak = selec- 
tion rule applies. In Ref. Il5|, Meier et al. observed pro- 
nounced decreases of the Raman shift in Silicon nanopar- 
ticles. A comparison with calculated values based on 



the phonon-confinement model showed the experimen- 
tally observed decrease of the Raman shift to be consider- 
ably stronger than predicted by the model calculations. 15 

Theoretically, the vibrational properties of nanostruc- 
tured silicon have also been studied. Heino has calculated 
phonon dispersion relations for silicon thin films using 
molecular-dynamics simulations^ whereas Saviot et al. 
have studied low frequency vibrational modes in silicon 
nanoparticles using an elastic continuum models Hu et 
al. and Valentin et al. employed empirical atomistic force 
models to calculate the VDOS of silicon nanoparticles by 
diagonalization of the dynamical-matrix i 19 i 20 The latter 
two studies did not include relaxation effects from the 
particle surfaces. 

In this work, we use molecular-dynamics simulations in 
combination with the modified embedded- Atom method 
(MEAM)2ii22 potential for silicon from Ref. [H] to cal- 
culate the vibrational density of states and selected zone- 
boundary phonon frequencies for silicon nanoparticles 
with diameters in the range 2.3 - 10.3 nm. These simula- 
tions include surface effects and the results show changes 
from the bulk density of states that are similar to those 
observed in metallic nanoparticlesi 2 ^— While the limita- 
tions of the potential do not allow us to make quantita- 
tive predictions, our calculations show clear trends in the 
developement of the particle VDOS as the particle size 
is reduced. In particular, we observe an increase of the 
VDOS at low frequencies, a shift of the high frequency 
modes to higher frequencies and a transfer of weight from 
the TO mode peak to intermediate frequencies. Our re- 
sults provide an explanation for the discrepancy observed 
in Ref. M- 



II. COMPUTATIONAL DETAILS 

The calculations of vibrational properties of Si 
nanoparticles presented in this work are all based on 
molecular-dynamics simulations using the MEAM poten- 
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TABLE I. Approximate diameters d of the model particles 
and the number of atoms N they contain. 



d(nm) 


2.3 


2.8 4.1 6.1 


8.2 10.3 


N 


329 


657 1803 5815 


14175 27943 



tial for Si from Ref. |23|. The advantage of the MEAM 
formalis m 21 ' 22 is that it is not only suitable for simula- 
tions of pure Si in the diamond structure but that it has 
also been used successfully in studies of composite sys- 
tems combining Si with other elements like, Ni,— Mb, 28 
or Aui22 A disadvantage of the potential is that it severely 
overestimates the phonon frequencies in silicon. 18 Qual- 
itatively the potential gives however a fair description 
of silicon and reflects the characteristic features of its 
phonon dispersion relations. Since in this work we are 
more interested in the qualitative changes of the VDOS 
than quantitatively exact frequencies and since we hope 
to exploit the potential's transferability in future work 
focusing on composite systems, we decided to use this 
potential despite its shortcomings. 

In order to study the influence of particle size on the vi- 
brational properties in Si nanoparticles, we constructed a 
set of model particles by cutting a spherical region out of 
an ideal diamond lattice. The diameters of the particles 
and the number of atoms they contain are summarized 
in Table |TJ In addition to the nanoparticle configurations 
we used a bulk-like configuration containing 1,000,000 Si 
atoms arranged on an ideal 50 x 50 x 50 cells diamond 
lattice. 

Prior to the calculation of vibrational frequencies, we 
carefully equilibrated our model systems at T = 300 K. 
In order to determine the equilibrium lattice constant of 
the bulk configuration at this temperature, the Parinello- 
Rahman technique was employed to perform simulations 
at constant pressure (p = 0). After the determination of 
the lattice constant, the configuration was re-equilibrated 
with a fixed volume. All simulations were carried out 
using the velocity Verlet algorithm with a time step At = 
0.5 fs. In simulations of the bulk configuration, periodic 
boundary conditions were applied. 

For the calculations of vibrational properties from 
our simulations we relied on the velocity-autocorrelation 
function. The total VDOS g{v) of a system of atoms 
is proportional to the Fourier-transform of the velocity- 
autocorrelation function averaged over all atoms. The 
normalized VDOS of an N atom system is thus given by 
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where v, denotes the velocity of the i-th atom. Infor- 
mation about individual vibrational frequencies can be 
obtained by this method if the velocities of the atoms 
are projected onto a plane wave. 

For a mode with wave- vector q and a normalized po- 
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FIG. 1. (Color online) Normalized vibrational density of 
states of crystalline bulk Si and nanoparticles of various sizes. 
The data have been broadened with a finite line-width of cor- 
responding to an energy of ~0.12THz 



larization vector p we calculated the function 
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where 
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is the projection of the velocities on the plane wave and 
r° denotes the average position of the i-th atom. For a 
crystalline system with translational symmetry the func- 
tion <?p should have sharp peaks at the positions of the 
phonons with wave- vector q that have a component along 
the polarization vector p. 

The vibrational modes of systems without transla- 
tional symmetry are not plane waves. For systems similar 
to perfect crystals (for example a finite piece cut out of 
an infinite crystal) one finds however that the vibrational 
modes are well described by phonon-like wave-packets. 
Such wave-packets appear in g^(f) as more or less broad- 
ened peaks whose maxima can still be used to determine 
the dominant frequency of the mode. The quality of the 
peaks in g q (i / ) — if there are clear, isolated peaks — can 
therefore be used to judge to what extent the phonon 
picture applies to the system in question. 



III. RESULTS 

In Fig. Q] we compare the VDOS that we obtained for 
three selected nanoparticles with the result of crystalline 
bulk Si. Since for the smaller particles the discrete na- 
ture of the spectrum becomes notable, we have broadened 
the VDOS data with a finite line-width corresponding 
to an energy of 0.5 meV (0.12 THz). As observed pre- 
viously in Ref. 0, the MEAM potential from Ref. [H 
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FIG. 2. (Color online) High frequency part of the normalized 
vibrational density of crystalline bulk Si and nanoparticles of 
various sizes. 
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FIG. 3. (Color online) Static structure factor along the line 
q = [qqq] for Si nanoparticles of various sizes. 



leads quantitatively to strongly exaggerated phonon fre- 
quencies. Qualitatively however, the results show a clear 
evolution of the VDOS as the particle size is diminished 
with the strongest changes occurring at the upper and 
lower end of the frequency spectrum. 

It can be seen from Fig.[T]that with decreasing particle 
size the large peak of transverse acoustic (TA) modes 
in the rage - lOThz is shifted to lower frequencies. 
The resulting increase of the VDOS is further amplified 
by the fact that the curves become more linear for the 
smaller particles. These two changes lead together to a 
significant increase of the low- frequency VDOS. There is 
however no notable loss of the height of the TA mode 
peak. 

In the high-frequency range the VDOS of the Si 
nanoparticles show a different behavior. As shown by 
Fig. [1] the strong peak of the TO modes around 25 THz 
loses its height rapidly with increasing particle size with- 
out a significant increase in peak width. In addition to 
this, the TO peak is slightly shifted to higher frequencies. 
This can also be seen in Fig. [2] where we show the high 
frequency flank of the TO peak without line broadening. 
This figure makes it clear that for the nanoparticle sys- 
tems this flank of the peak and the cutoff frequency are 
located at higher frequencies than in the bulk system. 
Due to the lower peak height the VDOS curves of the 
particles eventually cross the bulk line at frequencies be- 
tween 25 and 26.5 THz but it is clear from Fig. Q] that 
this is not due to a peak broadening but a shift to higher 
frequencies. The shift of the TO peak will be discussed 
in more detail below by looking at the frequencies of in- 
dividual vibrational modes. 

Interestingly we observe a non-monotonous behavior 
of the cut-off frequencies as shown by the fact that the 
cut-off frequencies of the Si329 and Si27943 particles are 
both lower than that of Siiso3- We interpret this as an in- 
dication that our smallest particles are below the limit of 
the scaling regime in which the physical properties scale 



monotonously with the particle size in a simple manner. 
This view is confirmed by other results discussed below. 

At intermediate frequencies, the VDOS data presented 
in Fig.Q]show only little differences between the nanopar- 
ticles and the bulk system. In particular, the shape and 
the position of the longitudinal acoustic (LA) modes peak 
at approximately 14 THz remains remarkably constant 
for all of our model systems. On the other hand, the 
position of the longitudinal optical (LO) mode peak at 
18 THz is slightly shifted to higher frequencies. In the 
regions to the left and to the right of the LO peak the 
VDOS increases with decreasing particle sizes. This in- 
crease accounts obviously for the loss of weight in the TO 
peak. 

Capillary pressure is a pressure in a material that is 
caused by curved surfaces or interfaces. The magnitude 
of this pressure depends upon the curvature of the surface 
and the materials surface or interface stress (which is in- 
dependent from the curvature). While this phenomenon 
is well known and easy to understand in the case of flu- 
ids, the situation is slightly more complex in the case of 
crystalline solid materials. A detailed overview of the ef- 
fects of curved surfaces in solid systems has been given 
by Wolf and Noziere^ 

In nanoparticles capillary pressure can have a notable 
effect on the lattice constant insdie the particle since the 
small size of the particles leads to extremely high surface 
curvatures. Such effects have previously been observed 
theoretically an experimentally in metallic nanoparticles 
(see for example Ref. [26| and references therein). In or- 
der to account for the effect of capillary pressure as well 
as surface relaxation effects on the lattice constant in 
our calculations of individual vibrational modes at high- 
symmetry points of the Brillouin zone, we have calculated 
the static structure factor of the model particles and ob- 
tained the lattice constant from the position of the (111) 
peak. 

Figure [3] shows the static structure factor along the 
[111] direction for all of our model particles. Except for 
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FIG. 4. (Color online) Frequencies of selected Si nanoparticle 
phonon modes as a function of the inverse particle radius. 
Data at 1/R = Onm -1 represent crystalline bulk values. For 
further discussions see main article text. 



the smallest particle, the maximum of the (111) peak 
shifts to larger values of q with decreasing particle size. 
This is the expected behavior since the smaller particles 
should have a larger capillary pressure (assuming a posi- 
tive surface stress for Si) and consequently a smaller lat- 
tice constant. The fact that the position of the maximum 
of the Si329 particle moves against the trend is another 
indication that this particle is outside the scaling regime 
for Si particles. It should be understood here that the 
lattice constants derived from g max are only meant to 
be used as approximate, average values since the crystal 
lattice in the nanoparticles is not homogenous. Capil- 
lary pressure compresses the lattice inside the particles 
whereas the surface stress has the inverse effect in the 
vicinity of the surface. 

With the average lattice constants derived from the 
static structure factor calculations, we calculated the fre- 
quencies of the vibrational modes corresponding to the T, 
X, and L point of the Brillouin-zone of the diamond lat- 
tice for all of our model nanoparticles using the method 
outlined in Sec. [TXJ The frequencies obtained from the 



maxima of the peaks in the projected density of states 
g^(v) are presented in Fig. @] Note that at the T point 
the LO and TO modes and at the X point the LA and 
LO modes are degenerate. For this reason the frequen- 
cies of these modes are denoted by ^LTo(r) and fLAo(X), 
respectively. 

For the larger particles and the bulk system the analy- 
sis procedure gave straightforward, unambiguous results 
indicating that the vibrational modes are sufficiently sim- 
ilar to plane waves. In case of the smaller particles Si329 
and Si657, however, the projected density of states some- 
times showed a group of similar peaks instead of a single 
peak. If in such a case it was still possible to identify one 
peak that was clearly most consistent with the rest of 
the data we marked the corresponding data point with 
an open symbol in Fig. 2] In the case of the ^lao(X) 
mode, however, it turned out to be impossible to iden- 
tify matching peaks for all particles smaller than Sissis 
so that these points were dropped from the figure. 

A special case is the z^ta(L) mode. For this mode, the 
projected density of states of all nanoparticles showed 
two clearly distinguishable peaks. This indicates that the 
degeneracy of this mode is lifted in the finite particles. 
For this reason we included two distinct data sets for this 
mode in Fig. @] 

The behavior of the individual frequencies is generally 
consistent with the changes in the total density of states 
shown in Fig. [I] In the bulk system, the TA modes at the 
L and X point define the upper and lower boundary of the 
flat plateau on top of the TA peak. The broadening and 
shift of this peak in the total VDOS compares well with 
the fact that ^ta(X) remains nearly constant, whereas 
the lower branch of the L-point TA mode shifts to lower 
frequencies for smaller particles. 

In the intermediate frequency regime, the sharp max- 
ima of the LA and LO mode peaks correspond to ^la(L) 
ad fLo(L) in the bulk case. Again, the very small shift of 
j/la(L) and the slightly stronger upward shift of i/lo(L) 
compare well with the changes in the total VDOS. This 
indicates that these modes are well described by the bulk- 
phonon picture. 

At high frequencies, ^to(X) and ^to(L) correspond 
to the lower and upper bound of the TO peak in the 
bulk case, whereas ^LTo(r) defines the cut-off frequency. 
The general increase of these frequencies with decreasing 
particle size (except for S1329) reflects the shift of the 
TO peak and the cut-off point in the same direction as 
shown in Fig. [1] and O Moreover, the data in Fig. @] 
support a slight sharpening of the TO peak; however the 
diminishing height of this peak in the VDOS makes it 
hard to verify this in Fig. [1] The slight decrease of the 
TO mode frequencies in the case of the smallest particle 
gives further evidence of the growing deviations from the 
bulk behavior in particles with diameters below 3 nm. 

Although the general trends of Fig. [4] agree well with 
the behavior of the VDOS, one has to be careful not to 
stretch these comparisons too far. As pointed out in the 
last paragraph, ^LTo(r) represents the cutoff frequency 
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FIG. 5. (Color online) Normalized local vibrational density 
of states in a Sissis nanoparticle compared to the density of 
states of bulk Si. 

in the bulk case. In case of Si32g, however, the calcu- 
lated value of ^LTo(r) = 26.80 THz is the same as in 
the bulk case, although Fig. [2] clearly shows that modes 
with higher frequencies are present in this cluster. This 
discrepancy is another hint that this particle is close to 
the limit of the regime where the phonon picture can be 
applied successfully. 

With decreasing size of a nanoparticles, the fraction 
of atoms that is located at or near its surface increases 
at the expense of the fraction of atoms in the interior of 
the particle. This means that the contribution of surface 
modes to the total VDOS becomes increasingly impor- 
tant for the smaller particles. 

In order to get insight into the influence of sur- 
face modes on the total VDOS, we have classified the 
atoms of one of our particles (Sissis) into three differ- 
ent groups: core atoms with a distance to the center 
d < r c = 2.35 nm, surface atoms with a coordination 
number Z < 4 and subsurface atoms with a perfect co- 
ordination Z — 4 but d > r c . The subsurface layer is 
necessary since atoms that are close to the surface gener- 
ally do not behave in the same way as atoms deep in the 
particle. The core, surface, and subsurface groups con- 
tain 2809, 835, and 2171 atoms, respectively. For each of 
the three groups we then calculated its contribution to 
the total VDOS of the particles. The resulting partial or 
local VDOS of the three regions in the Sissis particle are 
shown in Fig. [5] together with the VDOS of bulk Si. 

The partial VDOS shown in Fig. [S] reveal substan- 
tial differences between the three groups of atoms. The 
VDOS of the core atoms is very similar to the bulk curve. 
Below 15 Thz, in the range of the TA and LA modes, the 
two curves are nearly identical whereas the TA and TO 
peaks of the core atoms are slightly shifted to higher fre- 
quencies. This behavior of the core atoms agrees well 
with the trends of the bulk-like modes shown in Fig. 2] 

The VDOS of the surface atoms in Fig. [5] is markedly 
different from the bulk VDOS. At low frequencies a 
strong peak can be seen at frequencies below the TA peak 



of the bulk and core atoms. At intermediate frequencies 
the peak of the surface atoms LO modes appears to be 
shifted to lower frequencies and a broad band of modes is 
present around 20 Thz. At high frequencies, in the range 
of the TO modes of the bulk and core atoms, however, 
the surface VDOS is very low. 

A comparison of the trends observed in Fig. [T] with 
the surface VDOS in Fig. [5] shows that in those areas 
where the surface VDOS in Sissis is notably larger than 
the bulk VDOS, the total VDOS increases for the smaller 
particles. We therefore attribute the shift of the TA peak 
to lower frequencies, the increases of the VDOS at inter- 
mediate frequencies as well as the marked decrease of the 
TO peak to the growing influence of surface modes on the 
total VDOS. Not surprisingly, the subsurface atoms show 
a behavior somewhere between the two other groups. 



IV. SUMMARY AND CONCLUSIONS 

We have studied the vibrational properties of silicon 
nanoparticles with diameters in the range 2.3-10.3 nm us- 
ing molecular-dynamics simulations. Our results show 
clear changes of the VDOS as the particles become 
smaller. The strongest changes are visible in the low- and 
high-frequency part of the spectrum. At low frequencies, 
the VDOS increases as the particles become smaller and 
the TA mode peak is shifted to lower energies whereas 
at the upper end of the spectrum, the height of the TO 
mode peak is significantly reduced although the peak is 
shifted slightly (< 1 THz) to higher frequencies. Finally, 
an increase of the total VDOS is observed at interme- 
diate frequencies between the peaks of the LA and TO 
bands. 

Our calculations of the frequencies of individual vibra- 
tional modes and partial VDOS show that the observed 
VDOS changes are driven by two distinct mechanisms. 
The first mechanism is the growing fraction of surface 
atoms in the smaller particles which increases the contri- 
bution of surface modes to the total VDOS. As shown by 
the partial VDOS of the surface atoms in Fig.[5]this effect 
accounts for the shift of the LA peak to lower frequencies, 
the increase of the VDOS at intermediate frequencies and 
the notable decrease of the TO peak. 

While the surface modes lead to an overall downward 
shift of the VDOS, the picture changes if one looks at 
the frequency change of specific bulk-like phonon modes. 
As shown by Fig. 0] and confirmed by the behavior of the 
core atoms in Fig. [5] (the partial VDOS of the core atoms 
should not be affected by surface modes) , the frequencies 
of the bulk modes are either unchanged or increase as 
the particles become smaller. The most likely reason for 
this behavior is the presence of a capillary pressure inside 
the core of the particles that compresses and stiffens the 
lattice. 

In summary one can say that there are two different 
surface effects that affect the VDOS of Si nanoparticles 
in an opposing manner. While the appearance of surface 
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modes that are not present in bulk Si pushes the VDOS 
of the particles to lower frequencies, the frequencies of 
the remaining bulk modes are increased by the effect of 
capillary pressure. 

The increase of the VDOS at intermediate frequencies 
at the expense of the height of the TO peak that is caused 
by the surface modes provides a good explanation for the 
experimentally observed strong decrease of Raman shifts 
in silicon nanoparticles^ Calculations of Meier et alJ^ 
have shown that the measured Raman frequencies are 
lower than those predicted by the phonon confinement 
model. Their calculations, however, used the dispersion 
relations of bulk silicon and did not include any changes 
of the VDOS in the nanoparticles. Our results suggests 
that the transfer of weight from the TO peak to lower 
frequencies caused by the appearance of surface modes 
can account for the discrepancy between the experimen- 
tal results and the phonon confinement model. 

In general, our calculations agree well with previous 
studies of the vibrational properties of fee and bec metal 
nanoparticles <2i r 27 ' 31 ' 32 Similar to the observations in 
this study, the metallic nanoparticles show an increase 
and more linear behavior of the VDOS at low frequen- 
cies whereas the high frequency LA peak (note that there 
are no optical modes in the fee and bec lattices) is dimin- 
ished and shifted to higher frequencies. Partial density 
of states calculations have shown that the increase of the 
low-frequency VDOS and the loss of weight of the LA 
peak in copper nanoparticles can be attributed to modes 
localized at the particle surface s 26 ' 27 whereas the shift of 
the LA peak is caused by capillary pressure ! 25 ' 26 We con- 
clude that the TA and TO peaks in the Si nanoparticles 
undergo similar changes as the TA and LA peaks in the 
metallic systems. 



The structure of the calculated VDOS and the fact 
that we were able to derive phonon frequencies from the 
projected density of states shows that the phonon pic- 
ture can be successfully used to describe the vibrational 
modes in Si nanoparticles with diameters down to sizes 
below 4 nm. For particles with diameters of 3 nm and less 
our results indicate an onset of non-monotonous changes 
of physical properties as a function of the particle size. 
This means that the size of these particles is below the 
lower limit of the scaling regime where properties scale 
with the particle diameter in a simple manner. The vi- 
brational modes of particles in this range are still rather 
bulk-like, but deviations from the bulk picture become 
more pronounced. This finding agrees with Ref. [2(] where 
it was concluded that for silicon particles with diameters 
above 2.5 nm the VDOS becomes bulk-like. 

Unfortunately, the limitations of the MEAM potential 
employed in this work do not allow calculations of quan- 
titatively correct phonon frequencies. For this reason it 
would be interesting to verify our results with another 
model. While DFT calculations of the larger particles 
are currently out if reach, it might be possible to per- 
form such calculations for the smallest particles used in 
this study. 
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